# pip install plotly
import numpy as np
from matplotlib import pyplot as plt
import numpy.random as nprd
from scipy.stats import norm,cauchy
from IPython.display import clear_output
normal_d,cauchy_d = norm(),cauchy()
from time import time
import seaborn as sns
import plotly.graph_objects as go
import random as rd
A partir de l'inverse généralisée de la fonction de répartition
L'inverse généralisé de F est donné par : $F^-1 = -\dfrac{1}{\lambda} ln(1-u)$
def simulation_echantillons(n=1000,Lambda=2,norm=1):
tab_echantillons=[]
for i in range(n):
u=np.random.uniform()
tab_echantillons.append(((-1/Lambda)*np.log(1-u))/norm)
return tab_echantillons
plt.style.use('ggplot')
particules=simulation_echantillons(1000,2)
sum=np.sum(particules)
particules=[p/sum for p in particules]
plt.hist(particules,bins=100)
plt.title("Echantillons d'une variable aléatoire exponetielle")
Text(0.5, 1.0, "Echantillons d'une variable aléatoire exponetielle")
norm = lambda a : a/1000
x=np.linspace(0,4,1000)
f=2*np.exp(-2*x)
plt.figure()
plt.plot(x,f,label='densité de probabilité exponentielle pour lambda = 2')
particules=simulation_echantillons(1000,2)
plt.hist(particules,bins=100,density=True)
plt.ylim(0,2)
plt.legend()
<matplotlib.legend.Legend at 0x2d3e7aaba08>
On constate que les échantillons normalisés approchent l'aire sous la courbe de densité.
Pour N tend vers l'infini on montre que la somme des aires de l'histogramme tend vers l'aire sous la courbe.
Le problème de cette méthode est qu'elle demande une grande quantité d'achantillons pour estimer l'aire sous la courbe.
M=5
N_sample = 10000
def Accept_Reject(M=5,N_sample = 10000):
X_gen=[]
for j in range(N_sample):
u = nprd.uniform()
x = nprd.standard_cauchy()
while 1/M*normal_d.pdf(x)/cauchy_d.pdf(x) < u:
u=nprd.uniform()
x=nprd.standard_cauchy()
X_gen.append(x)
return(X_gen)
X_gen = Accept_Reject()
plt.style.use("ggplot")
plt.hist(X_gen,bins= N_sample//200)
plt.title(" Méthode Accept - Reject")
plt.show()
t=np.linspace(normal_d.ppf(0.0001),normal_d.ppf(0.9999),10000)
X=normal_d.pdf(t)
plt.plot(t,X)
plt.title(" Loi normale centrée réduite")
Text(0.5, 1.0, ' Loi normale centrée réduite')
plt.figure(figsize=[2*5.4, 1*2.8])
plt.hist(X_gen,bins= N_sample//200,color="grey",label="Sample Data",density=True)
plt.plot(t,X,color = "red",label="Densité de $\mathcal{N} (0,1)$",linewidth=2)
plt.title("Estimation de $X \sim \mathcal{N} (0,1)$ à l'aide d'une loi de Cauchy ")
plt.legend()
plt.show()
On constate qu'on approche bien une loi Normale centrée réduite. Ainsi, la méthode $\textit{Accept-reject}$ est une bonne méthode pour générer de nombres aléatoires d'une densité connue à l'aide d'un autre générateur de nombres aléatoires, aussi de densité connue.
On pourra discuter de l'influence de $M$ qui correspond à l'inverse de la probabilité que x soit accepté est égale.
plt.figure()
X_norm=normal_d.rvs(size=N_sample)
plt.hist(X_gen,bins= N_sample//100,color="blue",label="Normal With Cauchy",alpha = 0.4,density=True)
plt.hist(X_norm,bins= N_sample//100,color="red",label="Normal",alpha = 0.3,density=True)
plt.title("Estimation de $X \sim \mathcal{N} (0,1)$ à l'aide d'une loi de Cauchy et \n génération d'une loi normale centrée réduite")
plt.legend()
plt.show()
Qualitativement, on remarque bien que les deux histogrammes se supperposent et on en déduit qu'on a bien un générateur d'une loi normale à partir d'un générateur de nombres aléatoires suivant une loi de Cauchy
m = [1,2,4,8,16]
plt.figure()
T_exec=[]
for _,M in enumerate(m):
plt.figure(figsize=[2*5.4, 1.5*2.8])
t_deb=time()
X_gen = Accept_Reject(M)
t_exec = time()-t_deb
plt.hist(X_gen,bins= N_sample//20,color="grey",label="Sample Data",density=True)
T_exec.append(t_exec)
print(f"Temps d'exécution pour M = {M} vaut {t_exec} secondes" )
plt.plot(t,X,color = "red",label="Densité de $\mathcal{N} (0,1)$",linewidth=2)
plt.title(f" Méthode Accept - Reject avec $M={M}$")
plt.show()
plt.figure()
plt.plot(m,T_exec,color = "black")
plt.xlabel("$M$");plt.ylabel("Temps d'éxécution");plt.title("Temps d'éxécution en fonction de $M$")
plt.semilogx()
plt.show()
Temps d'exécution pour M = 1 vaut 1.6880464553833008 secondes
<Figure size 432x288 with 0 Axes>
Temps d'exécution pour M = 2 vaut 2.4730465412139893 secondes
Temps d'exécution pour M = 4 vaut 5.172055006027222 secondes
Temps d'exécution pour M = 8 vaut 10.372045755386353 secondes
Temps d'exécution pour M = 16 vaut 20.634044885635376 secondes
Plus M est grand, mieux on approxime la loi mais plus le temps d'exécution augmente
Dans l'exercice, au lieu de regarder les résulats à 1 dimension, on utilisera plotly pour regarder les résulats en 2 dimensions.
N_sample = 10000
def Box_Muller(N_sample = 10000):
X_gen,Y_gen=[],[]
for j in range(N_sample):
u1 = nprd.uniform()
u2=nprd.uniform()
R= -2*np.log(u1)
V=2*np.pi*u2
X_gen.append(np.sqrt(R)*np.cos(V))
Y_gen.append(np.sqrt(R)*np.sin(V))
return(X_gen,Y_gen)
X,Y = Box_Muller(10000)
plt.figure()
plt.scatter(X,Y)
plt.show()
plt.figure()
sns.kdeplot(x=X, y=Y,cmap="Blues",fill=True, thresh=0 )
plt.show()
Pour utiliser plotly, installez le à l'aide de $\textit{pip install plotly}$.
Si vous avez besoin de sauvegarder une image, utilisez fig.write_image("Image.png"), mais installez avant kaleido avec $\textit{pip install -U kaleido}$.
h,x1,y1,image = plt.hist2d(X,Y,bins=20,density=True)
fig = go.Figure(data=[go.Surface(z=h, x=x1, y=y1)])
# fig.write_image("Image.png")
fig.show()
plt.imshow(plt.imread("Image.png"))
<matplotlib.image.AxesImage at 0x2d3ecb0ca48>
On remarque qu'on génère bien une gaussienne en 2D à l'aide de la méthode de Box-Muller. Néanmoins, c'est une méthode plus couteuse en temps.
from scipy.stats import multivariate_normal
norm_2d = multivariate_normal([0,0],np.identity(2))
x,y=np.mgrid[-3:3:.1, -3:3:.1]
z= np.dstack((x, y))
fig = go.Figure(data=[go.Surface(z=norm_2d.pdf(z), x=x, y=y,opacity=0.7,colorscale="Greens"),go.Surface(z=h, x=x1, y=y1, opacity=0.4,colorscale='Reds',showscale=False)])
fig.write_image("Image3.png")
fig.show()
plt.imshow(plt.imread("Image2.png"))
<matplotlib.image.AxesImage at 0x2d3ec90a888>
plt.imshow(plt.imread("Image3.png"))
<matplotlib.image.AxesImage at 0x2d3ed0bad48>
La gaussienne 2D se confond avec l'histogramme 2D généré avec la méthode de Box-Muller. En relançant plusieurs fois le programme, on remarque qu'on conserve bien la coincidence entre l'histogramme et la surface de densité de la gaussienne bivariée.
Désormais, générons une loi normal de moyenne -3 et d'écart type 3.
def Box_Muller(N_sample = 10000,mu=-3,sigma=3):
X_gen,Y_gen=[],[]
for j in range(N_sample):
u1 = nprd.uniform()
u2=nprd.uniform()
R= -2*np.log(u1)
V=2*np.pi*u2
X_gen.append(sigma*np.sqrt(R)*np.cos(V)+mu)
Y_gen.append(sigma*np.sqrt(R)*np.sin(V)+mu)
return(X_gen,Y_gen)
X,Y = Box_Muller(N_sample = 5000,mu=-3,sigma=3)
from scipy.stats import multivariate_normal
norm_2d = multivariate_normal([-3,-3],3**2*np.identity(2))
x,y=np.mgrid[-10:4:.1, -10:4:.1]
z= np.dstack((x, y))
h,x2,y2,image = plt.hist2d(X,Y,bins=30,density=True)
fig = go.Figure(data=[go.Surface(z=norm_2d.pdf(z), x=x, y=y,opacity=0.7,colorscale="Greens"),
go.Surface(z=h, x=x2, y=y2, opacity=0.4,colorscale='Reds',showscale=False)])
fig.show()
On obtient ce résulat en remplaçant $\sqrt{R}\sin{V}$ par $\sigma\sqrt{R}\sin{V}+\mu$
from scipy.stats import multivariate_normal
def gen_N_mu_sigma(mu,sigma,d=1000):
n,X_gen=len(mu),[]
norm_nd = multivariate_normal(np.zeros(n),np.identity(n))
A = np.linalg.cholesky(sigma)
for j in range(d):
Z = norm_nd.rvs()
X_gen.append(Z@A+mu)
return(np.array(X_gen))
import seaborn as sns
X= gen_N_mu_sigma([10,1],[[0.4,0.2],[0.2,0.4]])
plt.figure(figsize=(10,10))
sns.kdeplot(x=X[:,0],y=X[:,1])
plt.show()
mu = np.array([0 ,50 ,100 , 50 ,100 ,200])
sigma = np.array([[11, 10, 5, 9, 4, 2]
,[10 ,13, 9, 15, 5, 3]
,[5, 9 ,15, 11, 3, 1 ]
,[9, 15 ,11 ,21, 6, 4]
,[4 ,5 ,3 ,6, 5, 1]
,[2 ,3 ,1 ,4, 1, 1]])
X= gen_N_mu_sigma(mu,sigma, 10000)
plt.figure(figsize=(10,5))
sns.kdeplot(X)
plt.show()
for i in range(6):
print(f"La variance de Z{i} est de {np.diag(sigma)[i]}")
La variance de Z0 est de 11 La variance de Z1 est de 13 La variance de Z2 est de 15 La variance de Z3 est de 21 La variance de Z4 est de 5 La variance de Z5 est de 1
Plus la variance est grande, plus la courbe de densité est étalée. C'est ce qu'on remarque sur le graphique ci dessus.
L'amplitude max est donnée par $$\dfrac{1}{\sqrt{2\pi}\sigma}$$
Ainsi, plus sigma est petit, plus l'amplitude max est grande.
import pandas as pd
sns.pairplot(pd.DataFrame(X,columns=[f"Z{i}" for i in range(6)]))
plt.show()
Indice, Corr = [],[]
for i in range(1,len(sigma)):
for j in range(i):
Corr.append(sigma[i,j]/(np.sqrt(sigma[i,i])*np.sqrt(sigma[j,j])))
Indice.append(f"{i},{j}")
pd.DataFrame(np.array([np.choose(np.flip(np.argsort(Corr)), np.array(Indice)),np.flip(np.sort(Corr))]).T,columns=["Indices","Corréalation"]).head(5)
| Indices | Corréalation | |
|---|---|---|
| 0 | 3,1 | 0.9078412990032038 |
| 1 | 5,3 | 0.8728715609439696 |
| 2 | 1,0 | 0.8362420100070909 |
| 3 | 5,1 | 0.8320502943378437 |
| 4 | 2,1 | 0.6445033866354896 |
On peut voir les différentes dispersions selon les diverses sous espaces.
Les variable $Z_3$ et $Z_1$ sont fortement corrélée et on peut le remarquer sur le graphique pairplot ci dessus. On a pu calculer les coefficients de corrélation : $$\rho_{i,j} = \dfrac{\Sigma_{i,j}}{\sigma{i}\sigma{j}}$$
def Bern(p,d=1000):
B_gen=[]
for j in range(d):
u = nprd.uniform()
B_gen.append(0 if u <p else 1)
return(np.array(B_gen))
B = Bern(0.7)
plt.hist(B,density=True)
plt.show()
nb_zeros = (len(B)-np.sum(B))/len(B)
print(f"La fréquence du nombre de zéros est de {nb_zeros}")
La fréquence du nombre de zéros est de 0.718
t_0 = 500
f = normal_d.pdf
c=0
x_t = [0]
while c <= 20000:
x,y = x_t[-1],nprd.uniform(-1,1)
rho = min(f(y)/f(x),1)
if rd.random() <= rho:
x_t.append(y)
else:
x_t.append(x)
c+=1
X = x_t[t_0:]
plt.style.use("ggplot")
plt.hist(X,density=True,bins=20)
t=np.linspace(normal_d.ppf(0.0001),normal_d.ppf(0.9999),10000)
X=normal_d.pdf(t)
plt.plot(t,X)
plt.title(" Loi normale centrée réduite")
Text(0.5, 1.0, ' Loi normale centrée réduite')
t_0 = 1000
f = normal_d.pdf
c=0
x_t = [0]
while c <= 20000:
x,y = x_t[-1],nprd.uniform(-20,20)
rho = min(f(y)/f(x),1)
if rd.random() <= rho:
x_t.append(y)
else:
x_t.append(x)
c+=1
X = x_t[t_0:]
plt.hist(X ,density=True,bins=30)
t=np.linspace(-5,5,10000)
X_f=normal_d.pdf(t)
plt.plot(t,X_f)
plt.title(" Loi normale centrée réduite")
Text(0.5, 1.0, ' Loi normale centrée réduite')
a1 = 1
a2 = 2
mu1 = 10
mu2 = -5
p = 0.3
import scipy.integrate as integrate
f= lambda x: p/(2*a1)*np.exp(-np.abs((x-mu1)/a1)) + (1-p)/(2*a1)*np.exp(-np.abs((x-mu2)/a2))
t_0 = 1000
c=0
x_t = [0]
while c <= 200000:
x,y = x_t[-1],nprd.normal((mu1+mu2-2)/2,3*(a1**2+a2**2))
rho = min(f(y)/f(x),1)
if rd.random() <= rho:
x_t.append(y)
else:
x_t.append(x)
c+=1
X = x_t[t_0:]
t=np.linspace(-20,20,10000)
X_f=f(t)
plt.hist(X ,bins = 1000,density=True)
a= integrate.quad(f,)
plt.title(" Loi de Laplace")
plt.plot(t,X_f)
[<matplotlib.lines.Line2D at 0x2d3f1749ac8>]
Prenons une variable aléatoire $X$ qui suit une loi normale de moyenne $\mu$ and d'écart type $\sigma$: $$X \sim \mathcal{N}(\mu,\,\sigma^{2})$$ de densité : $p(x) = \frac{1}{\sigma \sqrt {2\pi } }\exp({ - \dfrac{ (x - \mu )^2 } {2\sigma ^2 } })$$
et $$ \mu = \mathbb{E}(X) = \int_\mathbb{R} p(x) dx$$
Ainsi, on a un échantillage de Monte carlo avec $f(x) = x$ .
nb_particule = 100
X = nprd.normal(5,2,size=nb_particule)
mu_hat = np.mean(X)
print(f"La moyenne estimée vaut {mu_hat}")
La moyenne estimée vaut 5.130927466166048
from scipy.stats import norm
nb_particule = 100
V_M = []
for i in range (10000):
V_M.append( np.mean(nprd.normal(5,2,size=nb_particule)))
plt.hist(V_M,bins=200,density=True)
t=np.linspace(3,7,10000)
X_f = norm(5,2/np.sqrt(nb_particule)).pdf(t)
plt.plot(t,X_f)
plt.show()
La distribution de $V_M$ suit la distribution normal de moyenne $\mu = 5$ et de variance $\dfrac{4}{N}$
nb_particule = 1000
V_M = []
for i in range (10000):
V_M.append( np.mean(nprd.normal(5,2,size=nb_particule)))
plt.hist(V_M,bins=200,density=True)
t=np.linspace(3,7,10000)
X_f = norm(5,2/np.sqrt(nb_particule)).pdf(t)
plt.plot(t,X_f)
plt.show()
nb_particule = 10000
V_M = []
for i in range (10000):
V_M.append( np.mean(nprd.normal(5,2,size=nb_particule)))
plt.hist(V_M,bins=200,density=True)
t=np.linspace(3,7,10000)
X_f = norm(5,2/np.sqrt(nb_particule)).pdf(t)
plt.plot(t,X_f)
plt.show()
La variance $\dfrac{4}{N}$ diminue quand le nombre de particule augmente. Ainsi, l'estimation de $\mu$ s'affine quand on prend un plus grand nombre de particule.
import numpy as np
from scipy.stats import norm
from scipy.stats import gamma
import matplotlib.pyplot as plt
gamma.pdf(2, 2, 0, 1/3)
np.random.seed(123)
m0 = 5
s0 = 4
alpha = 3
beta = 0.5
nbChauffe = 20000
nbEchantillon = 2000 + nbChauffe
m = np.sqrt(s0) * norm.rvs(loc = m0, scale = 1, size = 1)
s = 1 / gamma.rvs(a = alpha, loc = 0, scale = 1 / beta, size = 1)
Y = norm.rvs(loc = m, scale = np.sqrt(s), size = nbEchantillon)
plt.hist(Y.flatten(),bins=1000)
plt.grid()
plt.show()
print('m : {}, s = {}'.format(m,s))
m : [8.61458824], s = [0.2877234]
from scipy.io import savemat
mesDonnees1 = {'Y' : Y,
'm0' : m0,
's0' : s0,
'alpha' : alpha,
'beta' : beta}
savemat('Gibbs.mat', mesDonnees1)
mesDonnees2 = {'m' : m,
's' : s}
savemat('GibbVraiesValeur.mat', mesDonnees2)
del m0
del s0
del alpha
del beta
del m
del s
del Y
from scipy.io import loadmat
matDonnees = loadmat('Gibbs')
matVraieDonnees = loadmat('GibbVraiesValeur')
Y = matDonnees['Y']
m0 = matDonnees['m0']
s0 = matDonnees['s0']
alpha = matDonnees['alpha']
beta = matDonnees['beta']
m = matVraieDonnees['m']
s = matVraieDonnees['s']
print(len(Y.flatten()))
22000
plt.figure()
plt.hist(Y.flatten(),bins = int(len(Y.flatten())/1000))
plt.grid()
plt.show()
from IPython.display import clear_output
X = np.zeros(shape = (2, nbEchantillon))
X[0, 0] = norm.rvs(loc = m0, scale = np.sqrt(s0), size = 1)
X[1, 0] = 1 / gamma.rvs(a = alpha, loc = 0, scale = 1 / beta, size = 1)
for i in range(1, nbEchantillon, 1):
M = (s0 / (nbEchantillon * s0 + X[1, i - 1])) * np.sum(Y) + m0 * X[1, i-1] / (nbEchantillon * s0 + X[1,i-1])
S = s0 * X[1, i-1] / (nbEchantillon * s0 + X[1, i-1])
X[0, i] = norm.rvs(loc = M, scale = np.sqrt(S), size = 1)
A = alpha + (nbEchantillon / 2)
u = (Y - X[0, i])
v = (Y - X[0, i]).reshape(-1, 1)
B = beta + 0.5 * (u @ v)
X[1, i] = 1 / gamma.rvs(a = A, loc = 0, scale = 1 / B, size = 1)
plt.subplot(2,1,1)
plt.plot(X[0, nbChauffe : nbEchantillon])
plt.grid()
plt.subplot(2,1,2)
plt.hist(X[0, nbChauffe : nbEchantillon])
plt.grid()
plt.show()
mHat = np.sum(X[0, nbChauffe : nbEchantillon]) / len(X[0, nbChauffe : nbEchantillon])
print('Estimation de m : ', mHat)
Estimation de m : 8.617519517345205
plt.subplot(2,1,1)
plt.plot(X[1, nbChauffe : nbEchantillon])
plt.grid()
plt.subplot(2,1,2)
plt.hist(X[1, nbChauffe : nbEchantillon])
plt.grid()
plt.show()
sHat = np.sum(X[1, nbChauffe : nbEchantillon]) / len(X[1, nbChauffe : nbEchantillon])
print('Estimation de s : ', sHat)
Estimation de s : 0.28727648846856063